## This script conducts SAOM models for "The Coevolution of Networks of Interstate Support, Interstate Threat and Civil War" 
## The starting point for the script was derived from Kinne & Bunte, "Guns or Money," BJPS 2018

## Since the estimation involves simulations, results may vary slightly from the published results.

library(reshape2)
library(network)
library(sna)
library(ggplot2)
library(gdata)
library(countrycode)
library(gdata)
library(parallel)
library(RSiena)


rm(list=ls(all=TRUE))

#dir <- "[Directory Information Here]"
dir <- "C:/Users/rapiduser/Box/Relation creation/JOPreplication/"


setwd(dir)

(n.clus <- detectCores() - 1)

## Import the network and monadic data in two separate objects
nets.total <- read.csv("relation_network.csv", header=TRUE, row.names=1)
mons.total <- read.csv("relation_nodal.csv", header=TRUE, row.names=1)
compchanges <- read.csv("compchanges.csv", header=TRUE, row.names=1)


## Set some parameters
yr <- sort(unique(mons.total$year))
yr.lo <- yr[-length(yr)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

## Make composition change list 
write.table(compchanges, "cc.dat",row.names=FALSE, col.names=FALSE,
            quote=FALSE)
changes <- sienaCompositionChangeFromFile("cc.dat")
names(changes) <- id

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variables into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly=FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Changing dyadic covariates, omitting the constant ones
dyad.effects <- d.cons[ !d.cons %in% c("contiguous") ]
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), varDyadCovar(array(do.call("c", get(x)), dim=c(N,N,nlo))))
}; rm(x)

## Constant dyadic covariates
for (x in d.cons[ d.cons %in% c("contiguous")]) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), varCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  support_mi_pot.net,
  threat_density_nbrs.net,
  intra_intens_nbrs.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net,
  
  changes

)

eff <- getEffects(netdata)

###Start with base model, and then add in network effects

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Estimate some models
#dir.create(paste(dir, "/Output", sep=""))
setwd(paste(dir, "/Output", sep=""))

tm <- Sys.time()
model <- sienaAlgorithmCreate(useStdInits=F, projname="Estimates.paper",
                              nsub=4, n3=5000, maxlike=F, firstg=.03,
                              modelType=c(support.net=1,threat.net=1))
print01Report(netdata, modelname="Report.base")

output <- siena07(model,data=netdata,effects=eff,batch=TRUE,
                  verbose=FALSE,useCluster=TRUE,nbrNodes=n.clus,returnDeps=FALSE)
save(output, file="output.paper_base")
xtable(output, file="Estimates_paper_base.tex", type="latex", digits=3)


keep(dir, sure=TRUE)


############################################
## Use above estimates to make forest plots
############################################

library(ggplot2)
library(grid)
library(gridExtra)


################Base model, Appendix Figure 1

setwd(paste(dir, "/Output", sep=""))
load("output.paper_base")

#dir.create(paste(dir, "/Figures", sep=""))
setwd(paste(dir, "/Figures", sep=""))

out <- cbind(output$theta,output$se)
rownames(out) <- output$effects$effectName

## First do support equation
out.support <- out[grep("support.net:",rownames(out)),]
for (i in 1:nrow(out.support)) { # Rescale means and SEs for forest plot with all SEs at 0.5
  ad <- 0.5/out.support[i,2]
  out.support[i,1] <- out.support[i,1] * ad
  out.support[i,2] <- out.support[i,2] * ad
  rm(ad)
}; rm(i)

#Rescale density and reciprocity
out.support[1,1] <- out.support[1,1] * 0.1
out.support[1,2] <- out.support[1,2] * 0.1
out.support[2,1] <- out.support[2,1] * 0.1
out.support[2,2] <- out.support[2,2] * 0.1

## Label rows and put everything in a data frame for ggplot
support.names <- c("Density (x 0.1)","Reciprocity (x 0.1)","Contiguity","Ideal Point Distance",
                   "Common IGO Memberships","j's Intrastate Conflict Intensity","i's Intrastate Conflict Intensity",
                   "j's Military Expenditures Ratio","i's Military Expenditures Ratio",
                   "j's Energy Consumption Per Capita","i's Energy Consumption Per Capita",
                   "j's Military Capabilities Change", "i's Military Capabilities Change",
                   "j's Democracy","i's Democracy","Joint Democracy")
dsupport <- data.frame(
  nms=support.names,
  x=out.support[,1],
  xlo=(out.support[,1] - (1.96 * out.support[,2])),
  xhi=(out.support[,1] + (1.96 * out.support[,2]))
)
p.cols <- vector()
for (ii in 1:nrow(dsupport)) {
  if (dsupport[ii,"xlo"] < 0 && dsupport[ii,"xhi"] > 0) {
    p.cols[ii] <- "#00b2a6" # Color for insignificant estimates
  } else {
    p.cols[ii] <- "#8519d7" # Color for significant estimates
  }
}; rm(ii)
dsupport$p.cols <- p.cols
#Adjust the order
dsupport <- dsupport[ c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16), ]

marg <- c(0.1,0.5,0.1,0)

p1 <- ggplot(dsupport) + geom_hline(yintercept=0, lty=2) +
  geom_linerange(mapping=aes(x=seq(1,nrow(dsupport)), ymin=xlo, ymax=xhi), color=dsupport$p.cols, linetype=1, size=1, alpha=0.5) +
  geom_point(mapping=aes(y=x, x=seq(1,nrow(dsupport))), size=3, shape=16, color=dsupport$p.cols) +
  scale_x_continuous(breaks=seq(1:nrow(dsupport)), labels=as.character(dsupport$nms)) +
  xlab("") + ylab("Rescaled estimates + 95% CIs") + ggtitle("Support Equation") + coord_flip() + 
  theme(
    panel.grid.minor.y=element_blank(),
    axis.text.y=element_text(color="black"),
    plot.margin=unit(marg,"cm"),
    plot.title = element_text(hjust = 0.5)
  )

pdf("SupportEstimates_paper_base.pdf",width=6,height=8,paper="special")
par(mar=c(0,0,0,0))
plot(p1)
dev.off()

keep(dir, out, output, sure=TRUE)

## Then do Threat equation
out.threat <- out[grep("threat.net:",rownames(out)),]
for (i in 1:nrow(out.threat)) { # Rescale
  ad <- 0.5 / out.threat[i,2]
  out.threat[i,1] <- out.threat[i,1] * ad
  out.threat[i,2] <- out.threat[i,2] * ad
  rm(ad)
}

#Rescale density and reciprocity
out.threat[1,1] <- out.threat[1,1] * 0.1
out.threat[1,2] <- out.threat[1,2] * 0.1
out.threat[2,1] <- out.threat[2,1] * 0.1
out.threat[2,2] <- out.threat[2,2] * 0.1

## Label rows and put everything in a data frame for ggplot
threat.names <-  c("Density (x 0.1)","Reciprocity (x 0.1)","Contiguity",
                   "Ideal Point Distance","Common IGO Memberships",
                   "j's Intrastate Conflict Intensity","i's Intrastate Conflict Intensity",
                   "j's Military Expenditures Ratio","i's Military Expenditures Ratio",
                   "j's Energy Consumption Per Capita","i's Energy Consumption Per Capita",
                   "j's Military Capabilities Change","i's Military Capabilities Change",
                   "j's Democracy","i's Democracy","Joint Democracy")
dthreat <- data.frame(
  nms=threat.names,
  x=out.threat[,1],
  xlo=(out.threat[,1] - (1.96 * out.threat[,2])),
  xhi=(out.threat[,1] + (1.96 * out.threat[,2]))
)
p.cols <- vector()
for (ii in 1:nrow(dthreat)) {
  if (dthreat[ii,"xlo"] < 0 && dthreat[ii,"xhi"] > 0) {
    p.cols[ii] <- "#00b2a6" # Color for insignificant estimates
  } else {
    p.cols[ii] <- "#8519d7" # Color for significant estimates
  }
}; rm(ii)
dthreat$p.cols <- p.cols
#Adjust the order
dthreat <- dthreat[ c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16), ]

marg <- c(0.1,0.5,0.1,0)

p2 <- ggplot(dthreat) + geom_hline(yintercept=0, lty=2) +
  geom_linerange(mapping=aes(x=seq(1,nrow(dthreat)), ymin=xlo, ymax=xhi), color=dthreat$p.cols, linetype=1, size=1, alpha=0.5) +
  geom_point(mapping=aes(y=x, x=seq(1,nrow(dthreat))), size=3, shape=16, color=dthreat$p.cols) +
  scale_x_continuous(breaks=seq(1:nrow(dthreat)), labels=as.character(dthreat$nms)) +
  xlab("") + ylab("Rescaled estimates + 95% CIs") + ggtitle("Threat Equation") + coord_flip() + 
  theme(
    panel.grid.minor.y=element_blank(),
    axis.text.y=element_text(color="black"),
    plot.margin=unit(marg,"cm"),
    plot.title = element_text(hjust = 0.5)
  )

pdf("ThreatEstimates_paper_base.pdf",width=6,height=8,paper="special")
par(mar=c(0,0,0,0))
plot(p2)
dev.off()

keep(dir, out, output, sure=TRUE)


##Finally, Intrastate Intensity Equation
rownames(out)[rownames(out)=="intra_intens.net linear shape"]<-"intra_intens.net: linear shape"
rownames(out)[rownames(out)=="intra_intens.net quadratic shape"]<-"intra_intens.net: quadratic shape"
out.intra <- out[grep("intra_intens.net:",rownames(out)),]
for (i in 1:nrow(out.intra)) { # Rescale
  ad <- 0.5 / out.intra[i,2]
  out.intra[i,1] <- out.intra[i,1] * ad
  out.intra[i,2] <- out.intra[i,2] * ad
  rm(ad)
}


## Label rows and put everything in a data frame for ggplot
intra.names <-  c("Linear Shape", "Quadratic Shape", "i's Military Expenditures Ratio",
                   "i's Energy Consumption Per Capita", "i's Military Capabilities Change", "i's Democracy")
dintra <- data.frame(
  nms=intra.names,
  x=out.intra[,1],
  xlo=(out.intra[,1] - (1.96 * out.intra[,2])),
  xhi=(out.intra[,1] + (1.96 * out.intra[,2]))
)
p.cols <- vector()
for (ii in 1:nrow(dintra)) {
  if (dintra[ii,"xlo"] < 0 && dintra[ii,"xhi"] > 0) {
    p.cols[ii] <- "#00b2a6" # Color for insignificant estimates
  } else {
    p.cols[ii] <- "#8519d7" # Color for significant estimates
  }
}; rm(ii)
dintra$p.cols <- p.cols
#Adjust the order
dintra <- dintra[ c(1,2,3,4,5,6), ]

marg <- c(0.1,0.5,0.1,0)

p3 <- ggplot(dintra) + geom_hline(yintercept=0, lty=2) +
  geom_linerange(mapping=aes(x=seq(1,nrow(dintra)), ymin=xlo, ymax=xhi), color=dintra$p.cols, linetype=1, size=1, alpha=0.5) +
  geom_point(mapping=aes(y=x, x=seq(1,nrow(dintra))), size=3, shape=16, color=dintra$p.cols) +
  scale_x_continuous(breaks=seq(1:nrow(dintra)), labels=as.character(dintra$nms)) +
  xlab("") + ylab("Rescaled estimates + 95% CIs") + ggtitle("Intrastate Conflict Intensity Equation") + coord_flip() + 
  theme(
    panel.grid.minor.y=element_blank(),
    axis.text.y=element_text(color="black"),
    plot.margin=unit(marg,"cm"),
    plot.title = element_text(hjust = 0.5)
  )

pdf("Intra_IntensEstimates_paper_base.pdf",width=6,height=4,paper="special")
par(mar=c(0,0,0,0))
plot(p3)
dev.off()

keep(dir, sure=TRUE)


